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ABSTRACT 

Accretion disks are three-dimensional, turbulent, often self-gravitating, magnetohydrodynamic flows, which 
can be modeled in detail with numerical simulations. In this paper, we present a new algorithm that is based 
on a spectral decomposition method to simulate such flows. Because of the high order of the method, we can 
solve the induction equation in terms of the magnetic potential and, therefore, ensure trivially that the magnetic 
fields in the numerical solution are divergence free. The spectral method also suffers minimally from numerical 
dissipation and allows for an easy implementation of models for sub-grid physics. Both properties make our 
method ideal for studying MHD turbulent flows such as those found in accretion disks around compact objects. 
We verify our algorithm with a series of standard tests and use it to show the development of MHD turbulnce 
in a simulation of an accretion disk. Finally, we study the evolution and saturation of the power spectrum of 
MHD turbulence driven by the magnetorotational instability. 

Subject headings: accretion disks — black hole physics — hydrodynamics — magnetohydrodynamics 
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1. INTRODUCTION 

Although the standard acc retion disk model was prop osed 
more than thirty years ago dShakura & Sunvaevlll973l) . the 
properties of turbulent angular momentum transport in ac- 
cretion disks are still not well understood. Shakura & Sun- 
yaev (1973) hypothesized in their original work that magnetic 
fields may be important in mediating the required angular mo- 
me ntum transport. However, it was not until the last decade 
that iBalbus & Hawleyl (Il991alfbh pointed out that the magne- 
torotational instability (MRI) generates turbulence and leads 
to transport of angular momentum in accretion disks. 

The non-linear evolution of the MRI a nd the generation 
of turbu lence was studied numerically by IBalbus & Hawleyl 
d!991bl) . following the earlier linear analysis of the insta- 
bility. Later, local numerical simulati ons were performed 
in the shearing box approximation (e.g..lHawlev et aDl!995b 
iBrandenburg et al.lll995UHawlev et alj|1996l) . aimed to study 
further the local properti es of three-dimen sional MRI, with 
and without stratification (IStone et al.lll996l) . The natural ex- 
tension of shearing box calculations, namely the cylindrical 
disks wi th vanish i ng ve rt ical gravitationa l force , were simu- 
lated by iHawlevi d2001l) ; lArmitage et~aT1 d2001l) to illustrate 
some important aspects of the turbulent transport, especially 
in the vicinity of the innermost stable circular orbit around 
a black hole. Finally, global numerical simulations of MHD 
disks have also been carried out for a variety of settings and 
physical conditions. 

All of the codes that have been used to study the properties 
of MRI-driven turbulence have been based on two types of 
differencing schemes. The first class of studies make use of 
the very successful scheme developed originally for the ZEUS 
code (e.g., Stone & Pringle 2001; Armitage et al. 2001; see 
Stone & Norman 1992a and 1992b for the ZEUS code) or 
schemes based on it (e.g., Hawley 2000; De Villiers & Haw- 
ley 2003; Steinacker & Papaloizou 2002; Igumenschev et al. 
2003). The other class of studies have used different conser- 
vative schemes (e.g., Koide et al. 1999; Gammie et al. 2003; 
Machida & Matsumoto 2003). Both types of finite differ- 
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ence methods allow for a stable and efficient implementation 
of solvers of the MHD equations. However, they also intro- 
duce a considerable amount of numerical dissipation to the 
problem. This is significant because most calculations have 
been performed for ideal MHD and, hence, it is this numeri- 
cal dissipation that allows for the MRI instability to saturate 
and the resulting turbulence to reach a dynamical steady state 
(see, however, Flemming, Stone, & Hawley 2000 where the 
assumption of ideal MHD is relaxed). As a result, the ki- 
netic and magnetic energies of different simulations saturate 
at different levels de pending on the resolution and the scheme 
dHawlev et al.ll 1999b . The effect of this shortcoming can be 
reduced, e.g., by increasing the resolution, by increasing the 
discretization order, or by using nu merical schemes that re- 
duce numerical diffusion (see, e.g.. iGardiner & Stonell2005l 
for an unsplit Godunov MHD code). 

In this third paper of the series, we address this issue by 
developing a version of our pseudo-spectral numerical algo- 
rithm to simulate three-dimensional MHD disks. Spectral al- 
gorithms are high order numerical methods, in which dynami- 
cal variables are evolved along orthogonal modes. For smooth 
functions, they require only ~ n grid points to accurately re- 
solve one wavelength, compared to ~ 16 grid points for finite 
difference schemes to reach the same accuracy. Moreover, the 
MHD equations can be evolved in time in terms of the vector 
potential A, preserving thus trivially the cylindrical character 
of the magnetic field. 

The high order of spectral methods makes them ideal 
for studying problems of magnetohydrodynamic turbulence, 
since they do not suffer from serious numerical dissipation. 
Moreover, spectral methods can easily incorporate models 
of sub-grid physics, such as the those involved in large - 
eddy simulation (see, e.g.. lBlackburn & Schm idt 2003]). The 
idea of large-eddy-simulations is to model approximately the 
small-scale structures, in stead of resol ving all features of a 
turbulent flow (see, e.g., ISagautll2004l) . with sub-grid mod- 
els that are based either on experiments or phenomenologi- 
cal descriptions of the small-scale physics. The large-eddy- 
simulation approach fits naturally to spectral methods be- 
cause the scale-separation operation is mathematically a fre- 
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quency low-pass filter, which is easily incorporated in a spec- 
tral scheme. 

In the next section we present the system of equations we 
solve. In particular, we describe our treatment of the induc- 
tion equation in terms of the vector potential, in which we in- 
troduce an appropriate gauge to avoid the unbounded (linear) 
growth in the vector potential in a Keplerian flow. In Sj3] we 
describe the numerical technique we use. In we present a 
series of tests to verify the implementation of our algorithm 
and conclude, in jj5j with an application of our method to 
the study of MHD turbulence in accretion disks driven by the 
magnetorotational instability. 

2. EQUATIONS AND ASSUMPTIONS 

Magnetohyrodynamics. — We consider three dimensional 
viscous, compressible, magnetohydrodynamic flows. The 
MHD equations contain four equations, namely, the continu- 
ity equation 

§+v-(H = o, (i) 

the momentum equation 

p— +p( v .V)v = -V(p+— J+ — (B-V)B+Vr+pg, (2) 
the energy equation 



hand side leads to a linear growth of A, in time. This growth 
will never saturate, because the mean of the product v^B z is 
always positive. Although the actual value of the magnetic 
field will not be affected, this growth will lead to a large round 
off error in A r if the MHD equations are integrated for a long 
time. This difficulty can be overcome by setting 



A 



'\l — = -2B z VGMr, 



(7) 



which suppresses the linear growth of the vector potential. In 
our algorithm, A is calculated dynamically from the values 
of B z (t,r) and v^,(f,r), where the over-bars indicate averages 
over the azimuthal and vertical directions of the disk, respec- 
tively. 

The analytical forms of the vari ous physical quan tities in 
equations (fl~|) — <[3j were given in IChanetal.ld200l . Here, 
we generalize them to three dimensions. The viscosity tensor 
(in Cartesian coordinates) is 



Tij 



2(/Z r + /i s )e,7+ /^r + Mb-^s (V ■ v)Sij 



where the strain-rate tensor etj is 



1 / dvi dvj 

2 \dxj dxi 



dE 



+ V ■ (E\) = -PV • v + <£>„ + — V ■ q — V ■ F, (3) The viscous dissipation rate is, therefore, 



and the induction equation 



OA c 2 

— =vx (VxA) + V 2 A + VA. 
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2(fr + p s )(eij) + /x r + yitb--At s (V-v 



(4) 



We denote by p the density, by v the velocity, and by E the 
thermal energy. In the momentum equation, P is the thermal 
pressure, r is the viscosity tensor, and g is the gravitational 
acceleration. In the energy equation, there are two dissipative 
terms, namely, the Ohmic dissipation and the viscous dis- 
sipation <£>„. We use q to denote the heat flux vector and F to 
denote the radiation flux. 

We write the induction equation in terms of the vector po- 
tential A, so that the magnetic field is given by B = V x A. 
The symbol a here represents the electrical conductivity and 
we define the microscopic resistance by r/ = c 2 /4-ira. The last 
term, VA, in the induction equation is a gauge source/sink 
term, the purpose of which we explain below. 

It is straightforward to show that equation (01) leads to the 
standard induction equation 



Finally, the Ohmic dissipation rate, $> v , is given by 

$„ = - = f |VxB| 2 . 

We again assume an ideal gas law so that 

3k B T 



E = p 
P = p 



2^wh ' 
p,m H ' 



(8) 
(9) 
(10) 

(11) 

(12) 
(13) 



— = Vx (vxB) + (V-r/V)B. 



(5) 



For the induction equation, we typically set the resistance r\ 
to zero so that the diffusion term in the induction equation 
vanishes. 

Gravity. — We solv e for the gr a vitatio nal acceleration, g, in 
a similar way as in IChan et alj d2006l) . We first define the 
gravitational potential ^ by 

g = -V*, (14) 

which is given by the volume integral 



Because V x (VA) = 0, VA does not affect the magnetic 
field. However, we retain this gauge term because, by proper 
choice, it can be used to suspend a non-physical (i.e., numeri- 
cal) linear growth in A and hence improve the accuracy of the 
scheme. To illustrate this, we consider a Keplerian disk with 
a constant vertical magnetic field so that the vector v x B has 
a non-zero r-component. By assuming r\ = 0, the potential 
form of the induction equation reduces to 



-g i^ld^ 



x-x 



(15) 



over all space. Rewriting equation (l5i in differential form, 
we obtain Poisson's equation 



V 2 * =47rGp, 



(16) 



dA r 

dt 



= B- 



GM dA 
r dr ' 



(6) 



where G and M are the gravitational constant and the mass 
of the central object, respectively. The first term on the right 



with satisfying the boundary condition 5 , (f,°o) = at all 
times. When simulating accretion flows, the computational 
domain is usually finite. Based on its linearity, we can 
decompose Poisson's equation into two parts, i.e., 



V 2 * int =47rG/) int , 



(17) 
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and 



V 2 *ext = 47rG/9 e 



(18) 



where pi nt denotes the mass density within the computational 
domain, which in our case is the disk density, and pext refers to 
external sources such as the central object and/or a companion 
star. The gravitational field is then given by 



g = gint + ge 



-V(* int +* ext ). 



(19) 



For the gravitational field of the ce ntral object, we use the 
pseudo-Newtonian approximation of iMukhopadhvavl (120021) 
for gext, which takes the form 

2 



gext 



£L f GM \ 
"73 {—) 



-2{a/c)^GMr/c 2 + (a/c) 2 
y/GMr/c 2 (r-r s ) + a/c 



(20) 

Here, rs = 2GM/c 2 is the Schwarzschild radius and a is a 
parameter related to the spin of the central object. 
In order to solve for self-gravity using equation ( fT7T > within 
we compute the integral 



*int(f, 



G r PmtC/y)^ (21) 

J®{3) x-x' 



as in iChan e t al. (2006). We then use equation ( fT9l ) to obtain 
the total gravitational field and use it in the momentum equa- 
tion. 

Subgrid Physics. — For homogeneous and isotropic turbu- 
lence, the ratio between the most energetic scale to the vis- 
cous length scale is proportional to , (Re i l A ), where Re is the 
Reynolds number. This scaling law suggests that the compu- 
tational cost is proportional to G(Re^ for three-dimensional, 
time-dependent simulations. Although spectral methods of- 
fer efficient ways to capture small scale features, they are still 
unable to resolve, even in the shearing box approximation, 
the expected dynamical range in accretion disks down to the 
molecular viscous length scale. 

In order to reduce the computational cost, we need to intro- 
duce an artificial cut-off to the probl em. Following the Large 
Eddy Simulation approach (LES, see Sagaut 2004^ for a very 
detailed introduction), we will use our numerical algorithm 
to solve the MHD equations for scales larger than this cut- 
off scale and introduce an approximate model for the smaller 
scales. Introducing such a model is required by the non-linear 
character of the MHD equations, which allow for coupling 
between the simulated large scales and the unresolved small 
scales. 

For any physical quantity /, we denote by / the filtered 
(i.e., resolved) function and define the sub-grid fluctuation by 



f=f~f. 



(22) 



We can then formally decompose any non-linear product of 
two physical quantities fg in physical space as 



where 



fg={f+f'){g+g') 
= fg + Lij + Cij + Rij , 



Lij=fg~fg, 



(23) 



Qj=fg'+f% 

«ij=TF- (24) 

If the quantities / and g are different components of the veloc- 
ity, these higher order correlations are the Leonard tensor, the 



cross-s tress tensor, a nd the Reynolds sub-grid tensor, respec- 
tively (Sagaut 2004). Note that the Reynolds sub-grid tensor 
is different from the (physical) Reynolds stress tensor. It ap- 
pears because of the presence of the artificial cut-off and is 
independent of the molecular viscosity. 

The basic idea of Large Eddy Simulations is to devise a 
model for the three tensors d24T i that captures the physics of 
sub-grid turbulence. This is beyond the scope of the current 
paper. Here, we will approximate the sum Ljj + Qj + Rjj in 
spectral space by the the spectral filter (Chan et al. 2005) 



CT/3 



exp 



llnel 



(25) 



where (3 is the order of the filter; e is the machine accu- 
racy, which is of order 10~ 15 for double precision floating 
point numbers; and n and N are the point index and num- 
ber of points in spectral space, respectively. Applying this 
filter once after every timestep is equivalent to adding a high- 
order (super) diffusion term in the dynamic equations. Al- 
though this approach is not based on any physical model, it 
has been known t o reproduce the large-scale proper t ies of tur- 
bulen t flows (e.g.. lKaramanos & Karniad akis 2000t IPasquettil 
l2005h . 

3. IMPLEMENTATION OF THE PSEUDO-SPECTRAL METHODS 

In our current implementation of the pseudo-spectral 
method in three-dimensions we use cylindrical coordinates 
because we are interested in the study of geometrically thin 
accretion disks. It is trivial, however, to alter the geometry 
of the domain of solution, when necessary, and use spherical- 
polar coordinates. 

Along the radial direction, we use a Chebyshev collocation 
method, whereas, for the azimuthal and vertical directions, 
we choose the Fourier basis. This implies periodic boundary 
conditions for both the azimuthal direction (which is natural) 
and the vertical direction (which needs to be justified for each 
specific application). The computational domain is 2#( 3 > = 
[r m m, '"max] x [—71", tt) x [ _ Z,Z) and we require r m i n > to avoid 
the coordinate singularity at the origin. 

Formally, we expand any physical quantity f(t,r, <p,z) as 

f(tM,z) = £f nml (t)T n (rW m ^ Mz/Z - 



(26) 



Here T n is the n-th order Chebyshev polynomial and r £ 
-l.ll is the standardized coordinate in the radial direction 



(IChan et aLll2005l) . Note that the frequently used Chebyshev- 
Gauss-Lobatto grid 



rk = cos 



for < k < N 



(27) 



has the property that 

F - h = Ov-i - r N « N~ 2 . (28) 
This gives the time-stepping constrain A? °< N~ 2 for 
hyperbolic (wave-like) equations and Af « N~ 4 for 
parabolic (diffusion-like) equations. To overcome this 
restriction, we use instea d the Kosloff-Tal-Ezer mapping 
dKosloff & Tal-Ezer l[T993h 



7" max 



arcsin(ar) 
arcsin(a) 



+ 1 



min 



arcsin(ar) 



i(a) 



(29) 



We choose the parameter a, which controls the regularity of 
the grid spacing, to be sech(|lne|/iV), where e is the ma- 
chine accuracy, to optimize the accuracy of the spatial deriva- 
tives dDon & Solomonoffll997l) . 
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We compute the numerical derivatives along the r- and in- 
directions as in IChan et all (I2005L 120061) . Hence, the radial 
derivative is given by the chain rule 

df = 1 df 
dr dr/dr dr ' 
where the derivative in the standardized coordinate is 

Jitup Mz/Z 



(30) 



n.m.l 



(31) 



We precompute analytically the derivative of the mapping, 
dr/dr. Here, we use /W to denote the Chebyshev coefficient 
of the radial derivative and compute it using the following 
three-term recursive relation 

fW n 
JNjnJ ~ U ' 



(32) 

&\ m ,l = 2NfN,m,l, (33) 

^!i/=£2, m ,/ + 2 (» +1 )/« + i,»,'' ( 34 ) 

where cq = 2 and c„ = 1 for n = 1,2, . . .,N. The azimuthal 
derivative is given by 



%=L imf nm ,(t)T n {f)e"^e 

t n,m.l 



un<p inlz/Z 



(35) 



Finally, the derivative along the z-direction is similar to the 
azimuthal derivative except for the extra normalization ir/Z, 
i.e., 

df ^ ml x 



dz 



I ^f m nl{t)T n (r)e'^e^l z . (36) 



n.m,l 



4. CODE VERIFICATION 



We have verified our numerical algorithm using a suit of 
test problems, some of which we present in this section. For 
a test particular to the three-dimensional hydrodynamics, w e 
adopt the free-falling dust ring test from IChan et alJ (120051) : 
for MHD , we study the mag netic braking of a rotating slab, 
following lStone et al.1 d 1 992b . 



4.1. An Advection Test: Free Fall of a Dust Ring 

Following IChan et all d2005t) . we use a free falling dust 
ring as an advection test of our three-dimensional algorithm. 
The computational domain is [0.2, 1 .8] x [— 7r, 7r] X [— 1, 1] with 
65 x 32 x 65 collocation points. The initial density is the 
Gaussian 



p = exp[-20(r-l) 2 -20z 2 ]. 



(37) 



In order to also test the advection along the z-direction, we set 
the initial velocity to v = (0,0, 1), which is equivalent with 
performing this test on a non-stationary Galilean frame. We 
assume that the gravitational acceleration is Newtonian and 
that of a central external object, i.e., that g = (-1/V 2 ,0,0). 
(Note that, for a cylindrically symmetric central object, the 
gravitational field should have been proportional to l/r.) We 
also neglect pressure and magnetic fields in this test. 

Th e analytical solutio n of this test is found in the same way 
as is IChan et al.l (120051) . Because the initial vertical velocity 
is non-zero, we have to replace z by z—v z t in order to capture 
the vertical motion. The analytical solution, therefore, reads 



E(f,r,0,z) = E (ro,z-v z f) — 



Ml 



+ r 




-1 o 1 

Horizontal axis x along <^ = 

FIG. 1 . — A comparison of the numerical to the analytical solution for the 
free-falling dust ring test problem discussed in the text. The gray-scale image 
depicts the density profile at different times in the simulation. The contour 
lines show the positions in the solution domain where the fractional error 
between the numerical and analytical solutions are 0.01% and 0.001%. In all 
plots, we also show the maximum error e raax for reference. 



where ro is found by solving implicitly the equation 



V2 



1 



f t = — sin I 2 arccos 
r o 2 



+ arccos , 



(39) 



In Figure Q] we plot the numerical density as gray-scale 
contours at = and for different times. The error between 
the numerical solution and the analytical solution is overplot- 
ted as a set of contour lines. The maximum error throughout 
the simulation is of order 10~ 4 . Note that our implementation 
of the free-streaming inner boundary condition does not intro- 
duce any significant errors caused by the artificial excitation 
or reflection of waves. 



(38) 
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FIG. 2. — A comparison of the numerical (circles) to the analytical (lines) solutions for the magnetic breaking test problem. The panels show a snapshot of the 
azimuthal velocity and magnetic field at time t = 50, for the discontinuous (DIC; upper panels) and the continuous (CIC; lower panel) initial conditions discussed 
in the text. For the discontinuous case, the numerical solution captures correctly the shock properties despite the Gibbs oscillations; for the continuous case, the 
numerical and analytical solutions are indistinguishable. 



4.2. A Test of Alfven Wave Propagation: Magnetic Braking of 
an Aligned Rotator 

The problem of magnetic breaking of an aligned rotator via 
the emission of nonli near, incompressible Alfven waves was 
solved analytically by Mouschovia s & PaleologovJ ( Il980h . It 
was then used by I Stone et all d 19921) to verify the ability of 
MHD algorithms to propagate transverse Alfven waves. In 
thi s test, we us e the sa me initial conditions as those described 
by lStone et al.l (119921) . namely, the discontinuous initial con- 
dition (DIC) and the continuous initial condition (CIC). Be- 
cause our algorithm is designed for compressible MHD flows, 
we discard the radial and vertical components of the momen- 
tum equation and set dP/dcj) = in the azimuthal component. 

We use the computational domain [0.2,1.8] x [— ir,w) x 
[-16, 16] with 33 x 32 x 512 collocation points. We set the 
initial density to 



P0 = 



1 

10 



kl>i 

otherwise 



(40) 



the initial magnetic field to B = (0,0, 1), and the initial veloc- 
ity to v = (0,rSlrj,0) where £Iq is the initial angular velocity. 
For the discontinuous case (DIC), we set the angular velocity 
to 

^° — | 1 , otherwise, 
whereas for the continuous case (CIC) we set it to 

, | Z |>1 

:j(l + cos7rz) , otherwise. 



O = 



(42) 



are 



given 



m 



The analytic solu t ions 

iMouschovias & Paleologould 19801) 

In Figured we compare the numerical to the analytical so- 
lutions for the two initial conditions at t = 50, right before the 
wave front passes z = 16. The numerical solution for the dis- 
continuous problem shows oscillations around the discontinu- 
ity (the Gibbs phenomenon), which are inherent to all spec- 
tral methods. However, the properties of the shock as well as 
those of the fluid around it are captured correctly in our nu- 
merical solutions. For the continuous problem, the numerical 
and analytical solutions are indistinguishable. 

5. TURBULENT MHD DISKS DRIVEN BY THE 
MAGNETO-ROTATIONAL INSTABILITY 

One of the advantages of spectral methods is the fact that 
they allow for an accurate control of numerical dissipation and 
hence t hey can be used to track accurately the stability of a 
flow. In lChanetal.1 ( 120051) we used our two-dimensional, hy- 
drodynamic spectral algorithm to successfully reproduce the 
Rayleigh sta bility criterion in a couette flow, even near the 
separatrix. In lChan et al.l ( 120061) we applied our spectral algo- 
rithm to self-gravitating disks and studied Toomre's criterion. 
In this section, we use our MHD spectral algorithm to study 
the properties of accreti on disks driven by the magnetorota- 
tional instability (MRI) dBalbus & Hawlevlll991al) . 

Any ionized and magnetized accretion disk is unstable to 
the MRI as long as the angular velocity of the flow is a de- 
creasing function of radius. If a cylindrical disk is threaded 
by a mean vertical magnetic field B z , all the waves along the 
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FIG. 3. — The power spectrum along the vertical direction of the mag- 
netic and kinetic turbulent energies in a Keplerian disk at r = 20, during the 
exponential growth of the MRI (at t = 200), for the simulation discussed in 
the text. The hash-filled area represents the limit of the numerical resolution. 
The vertical dotted line corresponds to the wavenumber of the most unstable 
mode as predicted by the analysis of the linear MRI. 

vertical direction with wavenumbers less than 



&mri — \/2q — , 

VA 

are unstable. Here, fl is the angular frequency, 



d\nr 



(43) 



(44) 



is a measure of the shear, and va = B z / y/4-np is the Alfven 
velocity for the background field. The growth rate of the in- 
sta bility depends strongly on the wavenumber and is given 
by (iBalbus & Hawleylll991al) 



7 : 



-(2-q)-k 2 +J(2-q) 2 +4k 2 



1/2 



The fastest growing mode has a wavenumber given by 
q /4 



^peak 



q 



VA 



(45) 



(46) 



We use our numerical algorithm to simulate the evolution 
of a Keplerian disk in a pseudo-Newtonian potential around 
a non-rotating object of mass M. For this calculation, we 
set G = c = 1 and all the distances and times in gravita- 
tional units. We solve the MHD equations in the computa- 
tional domain [3,43] x [-7r,7r) x [-2,2) using 257 x 64 x 32 
grid points. In order to justify the periodic boundary condi- 
tions along the vertical direction and avoid numerical prob- 
lems around the innermost stable circular orbit, we set the 
sound speed to w 0.2. We set the initial density to unity 
everywhere in the disk and thread the flow with a vertical 
magnetic field with a corresponding Alfven velocity equal to 
2 x 1(T 3 . In our dimensionless units, it takes a time of f w 60 
and t « 500, respectively, for the fluid elements to complete 
one orbit at r = 6 and r = 20. 

The initial exponential growth of the MRI offers another 
possibility to test the ability of our numerical method to cap- 
ture the properties of an MHD instability. For this reason, we 
discuss first the initial stages of the simulation and then the 
properties of the final state of saturated MHD turbulence. 

In Figure [3j we show the power spectrum of the magnetic 
and kinetic energies in the simulation, at a radius r = 20, dur- 
ing the initial linear regime of the instability. The similarity of 
the peak of the numerical power spectrum to the wavenumber 



FIG. 4. — The ratio of the r-<j> components of the Maxwell, and Reynolds 
stresses, as a function of radius in the accretion flow at time t = 2000. 
The error bars indicate the ~ 20% uncertainty in the value of each stress. 
Th e solid line shows t he result of the analytic calculation for the stress ra- 
tio IPessah et aU2006l) . 

of the most unstable mode of the linear MRI demonstrates 
that our numerical algorithm can reproduce, with an uncer- 
tainty comparable to the resolution, the wavenumber of the 
mode that shows the fastest growth rate (eq. 1146) ). 

The Maxwell and Reynolds stresses are amplified exponen- 
tially during the initial growth of the instability, with their ra- 
tios d etermined only by the value of local shear (Pessah et al. 
120061) . In Figure [4] we plot the ratio of the Maxwell to 
Reynolds stresses as a function of radius in the accretion 
di sk at time t = 200 an d compare it to the analytical solution 
of iPessah et al.l (120061) . The numerical and analytical ratios 
agree well, within the uncertainties of estimating appropriate 
averages. 

After a few orbital periods, turbulence is generated and the 
solution of the system of MHD equations becomes highly 
non-linear. Figure [5] shows the evolution of the density and 
magnetic energy in the accretion flow from the initial laminar 
state (first panel), through the time of the exponential growth 
of the MRI (second panel), to the final saturated state of MHD 
turbulence (last two panels). As also found in previous sim- 
ulations of MRI-driven turbulent accretion disks, the solution 
is highly variable, with large fluctuations in the various phys- 
ical quantities. Moreover, the magnetic and material stresses 
are not consistent with the prediction of the alpha model and 
are non-vanishing inside the innermost stable circular orbit. 
In this paper, we are particularly interested in the evolution 
and properties of the power spectrum of turbulence, which is 
a statistic that spectral methods are primarily suited to calcu- 
late. 

Figure [6] shows the evolution of the power spectrum along 
the vertical direction of the magnetic energy at cylindrical ra- 
dius r = 20. At early times, the magnetic energy increases 
exponentially because of the growth of the MRI, generating 
a spectrum of fluctuations that peaks around the wavenumber 
of maximum growth. At later times, interactions between dif- 
ferent modes lead to spreading of turbulent energy to nearby 
modes, generating a power-law spectrum of fluctuations. The 
index of the power spectrum of magnetic energy fluctuations 
is comparable to that of a Kolmogorov spectrum, even though 
the MHD turbulence is highly anisotropic. 

The power-law spectrum of magnetic energy fluctuations 
extends from the largest vertical scale of the simulation 
(which is the vertical extent of the domain of solution) to 
the smallest vertical scale (which is equal to the numerical 
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FIG. 5. — Four density (upper rows) and magnetic energy (lower rows) snapshots of the evolution of the MHD turbulence in an accretion disk driven by the 
magnetorotational instability. The face-on views show the two physical quantities at the midplane of the domain of solution, whereas the edge-on views are for 
an azimuth of <f> = 0. In each row, the first panel corresponds to a state very close to the initial laminar flow; the second panel shows the density during the 
exponential growth of the MRI; the last two panels show two different instants during the saturated state of the turbulence. 
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FIG. 6. — The evolution of the power spectrum along the vertical direc- 
tion of the magnetic energy at cylindrical radius r = 20, for the accretion 
disk simulation shown in Figure [5] The hash-filled area starts at the cut-off 
wavenumber of the spectral filter and represents the numerical resolution. 
The vertical dotted line corresponds to the wavenumber of the most unstable 
mode as predicted by the analysis of the linear MRI. 



resolution). Moreover, the integral of the power spectrum, 
which is equal to the total variance, depends on both scales. 
This is consistent with the fact that the saturation predictor 
found in numerical simulations of shearing boxes depends 
on bot h the vertical size of the box and the numerica l reso- 
lution (lHawley et al.lfl99l [l99l iPessah et al.ll2006bl) . It is, 
however, a result of two simplifying assumptions in our sim- 
ulation. The power spectrum extends to the largest vertical 
scale because we have neglected the vertical component of 
gravity and, therefore, the disk is not stratified. At the same 
time, the power spectrum extends to the smallest resolved 
scale because we have neglected Ohmic dissipation. Obtain- 
ing a realistic saturation predictor of MRI-driven turbulence 
in accretion disks will require numerical simulations of strat- 
ified flows with the largest possible dynamical range and an 
accurate model of sub-grid physics. 



C.-K. C. and D. P. acknowledge support from the NASA 
ATP grant NAG-5 13374. 



APPENDIX 

A. ADVECTIVE-CONSERVATIVE MIXED FORMALISM FOR MHD 

As in lChan et al.l d2005f) , we use the non-linear terms (v ■ V)v in their advective forms, whereas we use the terms that involve 
the density, p, and the energy, E, in conservative form. For the vector potential, we use an advective form in order to increase 
stability. In detail: 



9tp = _d^_d i (p^ 
r r 

a « Hi* \ a d r (rT rr -rP-rB 2 /87r) d^ r d z r zr P+B 2 /8n- 

O t V r = -V r a r V r - {OfhVr- v&) - V z O z V r + 1- 1- 1- 

r rp rp p rp 

BrdrBr B^B,-B^) B z 8 z B r 
~i : 1 : 1 : <~g r 
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4nrp 



4irp 
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4irp 
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B r d r B z B^d<j)B z B z d z B z 

4np 4irrp 4irp z ' 
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d,A r = Vj > B z -v z B c j, + r) 
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d,A z = v r B$- v^B r + 7) 



l -d r {rd r A r ) + ^dlA r + d}A r -^ 
d\A^ - *± + - r d r {rd r A^) + ^d 2 A 4 



-d r A 



r 



- d r (rd r A z ) + -=■ diA z + fif A, 
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(A6) 
(A7) 
(A8) 



The notation d t denotes partial deriv atives with respec t to tim e. We integrate this equations forward in time using a low storage, 
third-order Runge-Kutta scheme (see lChan et al.ll2005ll2006i for detail). 
We compute the magnetic field from the vector potential as 



B r = d J^-d z A, 
r 

B$ = d z A r - d r A z 

d r {rA<p) d^A r 

Dt = 



(A9) 
(A10) 

(All) 



The only non-trivial term here is d r (rA^), which is in conservative form. 
The viscosity tensor r i; in the above equation has the following general form 



Tij =2{p r + p s )e i j+ ( Pr + Pb-^s ) (V-v)£y. 



(A12) 



As we describe before, p r , pt,, and p s are the coefficients of radiative, bulk, and shearing viscosity. The strain rate tensor e,j 
written in cylindrical coordinate becomes 



e rr = d r v r 


(A13) 


ed>d> = + 
r r 


(A14) 


e zz = d z v z 


(A15) 


e r <p = = - 1 d r V0 - — + -a^v r J 


(A16) 


e<t,z = = 1 +o z v 1 


(A17) 


e-r = e rz = i (<9 z v r + d r v z ) . 


(A18) 



The viscous dissipation is 



= 2(fJL I + (JL t )(e i j) 2 + ( Air + A*b-^s ) (V ■ v) 2 



and the Ohmic dissipation rate is given by 



47T 



(A19) 



(A20) 



B. CODE PARALLELIZATION 



If the magnetohydrodynamic equations were linear, the parallelization of the spectral algorithm would be trivial, because the 
coefficient of each of the basis polynomials would evolve independently from the others. However, the non-linear terms in the 
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FIG. Al. — Schematic representation of the standard parallel FFT procedure in two dimensions (upper panel) and the alternate flip-flop procedure we have 
implemented in our algorithm (bottom panel). 



equations necessitate cross-processor communication. As discussed in lChan et al.l (2005), we incorporate the non-linear terms in 
the time-stepping by transforming the relevant physical quantities to coordinate space and evaluating the non-linear terms there, 
before transforming them back to spectral space. The efficient parallelization of our algorithm, therefore, relies on efficiently 
implementing a parallel version of the Fast Fourier Transform (FFT) algorithm for a multi-dimensi onal set of grid points. We use 
the Message Passing Interface (MPI) standard to parallelize our algorithm (see lGropp et al.|[l999alfbl) . 

In Figure [Al] we present schematically an example of our implementation of the FFT algorithm for transforming the values of 
a physical quantity on a two-dimensional grid along the x\ and %i directions. The standard algorithm for parallel Fast Fourier 
Transform starts with an one-dimensional ("slab") decomposition to distribute the data across different processors. We use solid 
lines in the figure to denote the distribution of grid points on different processors, which we label by PO, PI, P2, etc. It is clear 
from the upper-left panel of the figure that applying spectral methods along the X2 direction is straight forward. We can simply 
use the ordinary FFT for each row locally. However, for performing an FFT along the x\ direction, communication between the 
various processors is necessary. The standard method requires taking first a parallel transpose (PT), exchanging data between 
processors, and then a local transpose (LT), per processor, as illustrated in the upper panels of Figure [ATI These two processes 
exchange the directions of x\ and x-i and FFTs are then applied on each row separately. 

We developed an alternate ("flip-flop") procedure to increase the performance of parallel FFTs in mul ti-di mensions, by avoiding 
the time consuming step of obtaining the local transpose on each processor. The lower row of Figure [AT] illustrates the flipping 
procedure, which involves a block transpose (BT) and a parallel transpose (PT) operation. The flopping procedure follows the 
opposite routine. In these operations we can take advantage of the resulting memory layout and use the cache more efficiently. 
Because the order of the indices is preserved, the algorithm is easier to implement. More importantly, depending on the domain 
size and the number of processors we use, this approach can speed up the parallel FFT by ^20% in our 32-processor Beowulf 
cluster. 
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